#ifndef _EXCITATION_H
#define _EXCITATION_H

#include "GlobalDefines.h"
class Excitation
{
public:
	Excitation(){};
	virtual ~Excitation(){};
	virtual emx::Real3D GetIncident(emx::Real3D r, int n) = 0;
	bool SetTimeStep(emx::real_data delta)
	{
		Delta = delta;
		return true;
	};
	emx::real_data GetDelta(){return Delta;};
protected:
	emx::Real3D Rk, Pol;
	emx::real_data Delta;
};

class GaussianPulse:public Excitation
{
public:
	GaussianPulse(emx::real_data theta, emx::real_data phi, int pol, emx::real_data mag, emx::real_data pulsewidth, emx::real_data delay)
		:Mag(mag),PulseWidth(pulsewidth),Delay(delay)
	{
		emx::real_data Incident_Theta(theta*emx::PI/180), Incident_Phi(phi*emx::PI/180);
		emx::real_data TheSin = sin(Incident_Theta);
		emx::real_data TheCos = cos(Incident_Theta);
		emx::real_data PhiSin = sin(Incident_Phi);
		emx::real_data PhiCos = cos(Incident_Phi);
		//Incident unit vector
		Rk = -emx::Real3D(TheSin * PhiCos, TheSin * PhiSin, TheCos); 
		//electric field polarization
		if (pol == 0)
		{
			Pol = -emx::Real3D(TheCos * PhiCos, TheCos * PhiSin, -TheSin);
		}
		else if (pol == 1)
		{
			Pol = -emx::Real3D(-PhiSin, PhiCos, 0.0);
		}
	};
	virtual ~GaussianPulse(){};
	virtual emx::Real3D GetIncident(emx::Real3D r, int n)
	{
		double gamma = 4/PulseWidth * (emx::C*n*Delta - Delay - blitz::dot(r,Rk));
		return Pol*Mag*4 / (PulseWidth * sqrt(emx::PI)) * exp(-gamma * gamma);
	}
private:
	emx::real_data Mag, PulseWidth, Delay;
};

#endif